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Abstract — A  nonlinear  dynamical  signal  analysis  technique, 
recurrence  quantification  analysis  (RQA),  was  applied  to 
surface  electromyograms  (EMG)  recorded  during  a  series  of 
isometric  contractions.  None  of  the  ten  RQA  features  calculated 
adequately  related  the  EMG  to  the  force  level  so  principal 
components  analysis  was  applied  to  combine  these  features  into 
a  lower  number  of  variables.  Linear  regression  of  the  first 
principal  component  gave  similar  lines  for  each  subject. 
However,  the  error  was  too  great  for  these  lines  to  be  used  in 
predicting  force  from  the  principal  component. 

Keywords — EMG,  nonlinear  dynamics,  recurrence 
quantification  analysis,  principal  components  analysis. 

I.  Introduction 

A.  Aim 

The  aim  of  this  study  was  to  investigate  the  applicability 
of  a  nonlinear  dynamical  analysis  technique,  ‘recurrence 
quantification  analysis’  (RQA),  to  the  problem  of  examining 
the  electrical  activity  of  muscles  (electromyogram,  EMG) 
during  isometric  contractions.  Recently  [1],  it  was  observed 
that  features  obtained  through  RQA  may  vary  with  the  force 
produced  by  muscles.  It  has  also  been  reported  [2]  that  more 
meaningful  results  may  be  obtained  applying  principal 
components  analysis  (PCA)  to  a  set  of  RQA  features.  To 
investigate  this  for  EMG,  we  recorded  signals  from 
volunteers  during  muscular  contractions  at  specific  force 
levels  and  applied  a  joint  RQA-PCA  process. 

B.  Nonlinear  dynamics  and  embedding  theorem 

The  field  of  nonlinear  dynamics  is  concerned  with  the 
way  in  which  variables  describing  a  system  vary  with  respect 
to  each  other  instead  of  time.  If  the  system  is  described  by  d 
variables  xu  x2,  ...,  xj,  then  this  relationship  is  the  path 
traced  by  the  vector  [xi,  ...,  xj\ .  In  many  physical  and 
physiological  systems  however,  it  is  only  possible  to  measure 
a  single  variable.  For  example,  the  EMG  signal  might  be  the 
only  quantity  recorded  from  a  neuromuscular  system.  The 
dynamics  of  such  systems  can  be  modelled  by  applying  the 
‘embedding  theorem’  [3], 

The  embedding  theorem  is  based  on  the  concept  that  an 
observed  scalar  signal  is  a  ID  projection  of  the  dynamics  of 
the  system  which  are  represented  in  d  dimensions.  The 
dynamics  can  then  be  reconstructed  in  a  different  vector 
space  with  dE  dimensions,  using  only  the  data  in  the  scalar 
signal.  Suppose  that  the  recorded  data  is  a  sampled  and 
quantised  signal  s(n).  Then  embedding  vectors  y (n)  are 
constructed  from  successively  delayed  samples  of  the  signal: 

y(n)  =  \s(n),s(n  +  T), ...,«(«  +  ( dE  - 1)7' )] 


where  T  is  an  integer  multiple  of  the  sampling  period  and  dE 
is  the  chosen  dimension.  T  is  commonly  chosen  as  the  lag 
corresponding  to  the  first  zero  of  the  autocorrelation  function 
of  s(n),  although  it  has  also  been  argued  that  the  first  local 
minimum  of  the  auto  mutual  information  function  is  more 
appropriate  [4],  Little  difference  was  found  between  these 
two  measures  for  the  signals  recorded  in  this  study.  A  test  for 
false  nearest  neighbours  can  be  used  to  find  an  appropriate 
value  for  dE  [5].  A  vector  y F(n)  is  a  false  nearest  neighbour  to 
vector  y(n)  if  yF  is  the  closest  vector  to  y  when  they  are 
embedded  in  dimension  d,  but  much  further  away  in 
dimension  4+1.  Ideally,  dE  should  be  chosen  so  that  there 
are  no  false  nearest  neighbours;  in  practice,  <  1%  false 
nearest  neighbours  is  accepted  to  allow  for  the  effects  of 
random  noise. 

C.  Recurrence  plots  and  RQA 

Several  tools  have  been  developed  to  analyse  the 
dynamics  described  by  embedding  vectors  y («).  One  such 
technique  is  based  on  a  graphical  method  called  recurrence 
plot  analysis  [6].  A  recurrence  plot  shows  the  times  when  two 
vectors  y(i)  and  y (j)  are  close  to  each  other,  i.e.  when  y (j)  is 
within  a  distance  r  of  y(i).  If  y (j)  is  close  to  y(i)  then  point 
(i,j)  is  called  a  recurrent  point  and  a  plot  of  all  recurrent 
points  is  called  a  recurrence  plot.  Recurrent  points  form  a 
variety  of  patterns,  the  most  important  being  upwards 
diagonal  line  segments.  These  indicate  that  the  dynamics 
represented  by  a  series  of  vectors  are  later  repeated,  hence 
there  is  some  determinism  in  the  dynamics.  An  example  of  a 
recurrence  plot  is  shown  in  Fig.  1.  The  plot  is  symmetrical 
about  i  =  j  because  if  (i,j)  is  a  recurrent  point  then  from  the 
above  definition,  (j,  i )  is  also  a  recurrent  point.  Although  the 
way  in  which  recurrent  points  are  determined  seems 
straightforward,  it  is  difficult  to  obtain  much  useful 
information  from  the  plot  by  visual  inspection 

Recurrence  plots  are  most  usefully  described  using  a  set  of 
features  collectively  known  as  recurrence  quantification 
analysis  (RQA).  These  are  [6,  7]: 

•  REC  (%  recurrence),  the  percentage  of  the  recurrence 
plot  covered  by  recurrent  points; 

•  DET  (%  determinism),  the  percentage  of  recurrent  points 
contained  in  upwards  diagonal  line  segments; 

•  D/R,  the  ratio  of  DET  to  REC; 

•  Lmax,  the  maximum  length  of  upwards  diagonal  line 
segments; 

•  ER  (entropy  of  recurrence),  the  entropy  of  the 
distribution  of  lengths  of  upwards  diagonal  line 
segments; 

•  TLR  (trend  of  local  recurrence),  a  measure  of  the  change 
in  density  of  recurrent  points  away  from  the  line  i  =  j; 
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Fig.  1.  Recurrence  plot  of  one  second  of  a  surface  EMG  signal  recorded 
from  subject  2  at  30%  MVC  force  level. 

•  Ti  (mean  recurrence  time  of  the  first  type),  the  average 
time  taken  for  a  point  in  the  embedding  space  to  be 
‘revisited’; 

•  T2  (mean  recurrence  time  of  the  second  type),  the 
average  time  taken  for  a  point  in  the  embedding  space  to 
be  ‘revisited’  excluding  times  of  one  unit. 

In  addition  to  Lmax,  we  also  calculated  the  mean  and  median 
lengths  of  upwards  diagonal  line  segments,  Lmean  and  Lmed. 

The  TLR  measure  can  be  useful  when  applying  other  signal 
processing  methods  that  demand  stationarity,  since  stationary 
signals  have  a  TLR  value  close  to  zero. 

II.  Experimental  Methods 
A.  Data  Acquisition 

Eight  healthy  male  subjects  aged  18-36  years  (mean  = 
24)  volunteered  for  this  study.  Skinfold  thickness  as 
measured  over  the  right  biceps  brachii  muscle  was  2.5-6  mm 
(mean  3.5  mm).  Subjects  were  seated  and  instructed  to  flex 
their  right  wrist  against  a  cantilever  beam  with  their  upper 
arm  vertical,  elbow  flexed  at  90°  and  forearm  supine  (palm 
up).  The  beam  was  instrumented  with  four  strain  gauges  in 
full  bridge  arrangement.  Surface  EMG  was  recorded  from  the 
right  biceps  brachii  muscle  using  Noro-Trode™  adhesive  dual 
Ag-AgCl  electrodes  (Myotronics-Noromed,  Inc.)  and  63  dB, 
custom-built  10-500  Hz  EMG  amplifiers.  Prior  to  attaching 
the  electrodes,  the  recording  site  was  prepared  by  wiping 
with  an  isopropyl  alcohol  swab  to  dissolve  skin  oils,  shaving 
hair  from  the  recording  site  and  rubbing  the  skin  20  times 
with  800-grade  silicon  carbide  paper  to  thin  the  keratin  layer 
and  thereby  decrease  the  skin’s  resistance.  Data  was  recorded 
digitally  at  1000  samples  /  s  using  a  12-bit  data  acquisition 
card  with  input  limits  ±  5  V  (National  Instruments  PC-LPM- 
16PnP). 

Subjects  viewed  the  output  from  the  strain  gauge  bridge 
along  with  a  target  line  on  a  computer  monitor.  The  maximal 


voluntary  contraction  (MVC),  taken  as  the  maximum  force 
achievable  at  the  wrist,  was  determined  by  instructing  the 
subjects  to  follow  an  increasing  target  as  far  as  possible. 
EMG  was  then  recorded  while  the  subjects  attempted  to 
maintain  0,  10,  30,  50  and  70%  of  the  MVC  over  7  s. 
Subjects  rested  for  5-7  minutes  between  trials  in  order  to 
avoid  fatigue  effects. 

B.  Analysis 

RQA  was  performed  on  an  epoch  of  1  s  (1000  points) 
selected  from  each  recording.  Embedding  parameters  T,  cIe 
and  the  distance  threshold  r  were  different  for  each  subject, 
but  the  same  parameters  (calculated  from  the  70%  MVC 
recordings)  were  used  for  all  recordings  from  a  given  subject. 
T  was  taken  as  the  lag  corresponding  to  the  first  minimum  of 
the  auto  mutual  information  function,  calculated  as  described 
in  [8].  dE  was  chosen  as  the  lowest  dimension  for  which  the 
percentage  of  false  nearest  neighbours  was  <  1%  [4].  T  was  4 
samples  for  each  subject  and  values  for  dE  were  15-23  (mean 
=  17).  To  calculate  threshold  r,  first  values  corresponding  to 
30%  of  the  maximum  distance  between  embedding  vectors  at 
each  force  level  were  calculated,  r  was  then  chosen  halfway 
between  the  maximum  and  minimum  values. 

Since  the  ten  RQA  features  have  different  units,  each  was 
divided  by  its  standard  deviation  before  applying  PCA.  This 
‘standardised’  data  was  then  projected  onto  the  principal 
components  for  further  analysis. 

III.  Results 

A.  Individual  RQA  features 

Plots  of  the  changes  in  RQA  features  vs  force  levels  were 
not  encouraging.  Overall,  the  three  line  lengths  Lmax,  Lmean 
and  Lmed  decreased  with  increasing  force,  as  did  REC.  The 
ratio  D/R  increased,  but  there  was  little  overall  pattern 
otherwise.  Fig.  2  shows  plots  of  REC,  D/R  and  ER.  REC  and 
D/R  possibly  follow  power  laws,  but  it  is  difficult  to  tell  for 
certain  given  the  low  number  of  force  levels.  Certainly  there 
is  little  to  suggest  there  is  a  more  useful  relationship  between 
these  features  and  force  than  there  is  between  RMS 
amplitude  of  the  EMG  and  force. 

B.  RQA  features  combined  via  PCA 

It  was  found  that  the  total  variability  of  the  standardised 
features  was  completely  explained  by  four  principal 
components,  with  the  first  two  components  explaining  79- 
96%  (mean  86%).  Generally,  the  three  largest  weights  for  the 
first  component  (PCI)  corresponded  to  REC,  Lmean  and  Lmed; 
the  three  largest  weights  for  the  second  component  (PC2) 
corresponded  to  D/R,  ER  and  T2  or  TLR. 

Fig.  3  shows  the  RQA  features  data  for  all  subjects 
projected  onto  PCI  and  PC2  (i.e.  the  linear  combination  of 
the  RQA  features  indicated  by  the  weights  in  PCI  and  PC2), 
along  with  the  mean  results.  PCI  generally  decreases  with 
increasing  force,  although  not  at  all  uniformly.  PC2  generally 
has  a  local  minimum  at  30%  MVC,  except  for  one  subject 
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Force  (%MVC) 


Fig.  2.  Variation  of  some  RQA  features  with  force  (for  all  subjects). 

with  a  minimum  at  50%  MVC.  More  striking  is  the 
projection  of  the  data  onto  the  PC1-PC2  plane  (Fig.  4)  which 
shows  that  as  force  increases,  a  U-shaped  trajectory  is  traced 
out  clockwise  in  the  plane. 

Another  surprising  result  was  obtained  by  applying  a 
least-squares  linear  regression  to  the  projection  of  data  on 
PCI  (Fig.  5),  with  the  slopes  of  the  lines  for  each  subject 
appearing  similar  (between  -0.061  and  -0.088  units/%MVC, 
mean  =  -0.076  units/%MVC)  and  the  lines  intersecting  at  the 
32%  MVC  force  level  where  PCI  is  0. 

IV.  Discussion  and  Conclusion 

An  unfulfilled  goal  of  EMG  signal  processing  is  to  find  a 
definitive  relationship  between  the  EMG  and  force  produced 
by  the  muscle,  or  at  least  the  torque  about  a  joint  since  this  is 
easier  to  measure  than  muscular  force.  Such  a  relationship 
would  make  it  possible  to  determine  the  force  or  torque  from 
the  EMG.  While  some  of  these  results  are  interesting, 
particularly  the  behaviour  in  the  PC1-PC2  plane  and  the 
closeness  of  the  slopes  of  the  regression  lines  for  projections 
on  PCI,  it  is  unlikely  that  they  are  particularly  useful  in  terms 
of  this  goal.  The  difference  between  the  projected  RQA  data 
values  and  the  values  along  the  lines  of  best  fit  (Fig.  5)  are 
too  great  for  the  lines  to  be  used  in  predicting  force  levels. 
Also,  the  projections  on  PCI  were  only  monotonically 
decreasing  for  3  out  of  the  8  subjects  (Fig.  3).  Due  to  the  low 
sampling  of  this  curve  (only  five  force  values),  it  is  unclear 
whether  or  not  this  is  due  to  errors  for  the  other  subjects. 


PCI  PC2 


0  20  40  60  0  20  40  60 
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Fig.  3.  RQA  feature  data  in  terms  of  the  first  two  principal  components 
(PCI  and  PC2).  Dashed  lines:  results  for  individual  subjects.  Solid  lines: 
mean  result. 


PCI 


Fig.  4.  RQA  feature  data  projected  onto  the  PC1-PC2  plane.  Dashed  lines: 
results  for  individual  subjects.  Solid  lines:  mean  result. 

The  low  sampling  of  force  values  also  makes  the  task  of 
performing  PCA  an  ill-posed  problem,  since  there  are  more 
RQA  features  than  force  levels.  So  differences  between 
projections  on  the  components  for  different  subjects  may  be 
reduced  if  more  force  levels  were  examined.  Of  course,  this 
would  require  subjects  to  volunteer  much  more  of  their  time. 

Another  limitation  in  attempting  to  use  these  results  is 
that  this  analysis  was  performed  with  static  forces,  over  1  s 
epochs  during  which  the  EMG  signal  was  roughly  stationary 
(as  indicated  by  ITLRI  <  0.008  for  all  signals).  If  the  force 
changes  during  the  analysed  epochs  then  it  becomes  even 
more  difficult  to  relate  RQA  variables  or  a  linear 
combination  of  them  to  specific  force  levels.  But  in  another 
study  where  RQA  was  used  on  EMG  signals  for  studying 
fatigue  [9],  the  feature  DET  changed  earlier  than  other 
parameters  after  a  step  load  increase.  Whether  the  magnitude 
of  the  change  in  DET  was  necessarily  related  to  the 
magnitude  of  the  change  in  the  load  was  not  reported,  nor 
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Fig.  5.  Linear  regression  of  RQA  feature  data  projected  onto  PCI.  The 
actual  values  for  each  subject  are  also  shown  (cf.  fig.  3). 

were  the  effects  of  different  loads.  Interestingly,  and  in 
contrast  with  the  results  shown  in  Fig.  2,  there  was  little 
change  in  REC  in  that  study. 

One  problem  that  has  so  far  been  ignored  is  that  only  the 
EMG  of  m.  biceps  brachii  has  been  considered  as 
contributing  to  the  torque  about  the  elbow  joint,  when  in  fact 
there  are  at  least  two  other  main  elbow  flexor  muscles  that 
also  contribute.  It  is  unlikely  that  EMG  will  be  successfully 
related  to  the  torque  about  a  joint  until  methods  are 
developed  that  incorporate  the  activities  of  all  the  muscles 
contributing  to  that  torque. 
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